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I Abstract. The green alga Chlamydomonas swims with synchronized beating of its 

two flagella, and is experimentally observed to exhibit run-and-tumble behaviour 

'.^.j' similar to bacteria. Recently we studied a simple hydrodynamic three-sphere model 

^ of Chlamydom,onas with a phase dependent driving force which can produce run-and- 

Z/1 tumble behaviour when intrinsic noise is added, due to the non-linear mechanics of the 



cd 



system. Here, we consider the noiseless case and explore numerically the parameter 

space in the driving force profiles, which determine whether or not the synchronized 

I state evolves from a given initial condition, as well as the stability of the synchronized 

'^ state. We find that phase dependent forcing, or a beat pattern, is necessary for stable 

S synchronization in the geometry we work with. 
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PACS numbers: 05.45.45.Xt, 87.16.Qp, 47.63.-b 
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1. Introduction 

Microorganisms swim in the low Reynolds number regime where viscous forces dominate, 
inertia is negligible and the familiar propulsion methods of larger organisms become 
ineffective fllfl- Fluid flow is governed by the Stokes equation, which is time reversible. 
A necessary condition on a periodic swimming stroke in order to achieve net propulsion 
is that it is non-time reversible [5]. Inspired by sperm cells, which achieve propulsion 
by propagation of bending waves through their fiagellum, Taylor demonstrated that 
propulsion is possible in a viscous environment by studying the propagation of waves on 
an infinite sheet IgI. Purcell showed that a swimmer needs at least two compact degrees 
of freedom to break the time reversal symmetry and achieve net propulsion [I] . 

Many microorganisms swim using fiagella (^ ; there are two fundamentally different 
types of fiagella: bacterial fiagella and eukaryotic fiagella (or cilia). Eukaryotic fiagella 
form bends when microtubules on one side of the fiagella 'walk' or 'slide' along the 
microtubules on the other side (sl. The propagation of bends allows the fiagella to form 
beat patterns that can break the time reversal symmetry. For example, the first half 
of an individual cilium's beat cycle, called the power stroke, has the cilium sticking out 
and pushing the fluid, while the second half, called the recovery stroke, has the cilium 
bent as it returns to its original position jol. 

Our understanding of propulsion at low Reynolds number has been developed by 
theoretical model microswimmers. Lighthill demonstrated a model that can achieve 
net propulsion by studying periodic shape deformations of a nearly spherical swimmer, 
showing that the swimming velocity is at most of the order of the square of the amplitude 



of the deformations 10 . Purcell's three-link swimmer was studied by Becker et al., 
who determined the swimming direction and velocity for different angle amplitudes and 
relative link lengths fl 11 . A useful one-dimensional model is the linear three-sphere 



swimmer, where three beads are connected by two rods that change length with a non- 



reciprocal pattern 12 . Dreyfus et al. studied a rotational analogue of the three-sphere 



swimmer [13]. Avron et al. presented a more efficient swimmer consisting of a pair 



of bladders which exchange their volume and vary the distance between them 14 
There have been several experimental realisations of artificial low Reynolds number 
swimmers 



15 17 



When two sperms swim close to each other, their tails beat in synchrony 18 
and Taylor studied this using his waving sheet model with hydrodynamic interactions 
[gI. Coordinated beating of fiagella or cilia is important for a range of processes 



including motility, efficient pumping of fiuid and symmetry breaking in developing 



embryos 18-22 . Theoretical and experimental models have been studied to show that 



synchronization can occur through hydrodynamic interaction and that it is relevant 



to bacterial swimming and pumping by arrays 23 44 . Flagellar synchronization is 



observed in Chlamydomonas, a unicellular green alga that swims using two fiagella that 



beat with a breaststroke pattern 45-47 . The cell has diameter ~ 10/im and swims 
with velocity lOO/im/s so the Reynolds number is Re ~ 10~^ and inertia is negligible. 
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Figure 1. A three sphere model of Chlamydomonas. The left and right beads 
represent the flagella and move on circular trajectories in the cell frame. The back 
bead represents the cell body. The scaffold is shown to define the reference frame of 
the cell which has origin Rq in a lab frame, but it does not interact with the fluid. 
The green underlay is a schematic of a Chlamydomonas cell. 



During normal swimming, the flagella beat in synchrony. These periods of synchrony 
are interrupted by periods of asynchronous beating and during these asynchronies, there 
is a large change in the cell's orientation 48 , 49 . This is analogous to run-and-tumble 
behaviour observed in bacteria. 



Simple models have helped us understand better the intricacies of low Reynolds 
number swimming and hydrodynamic synchronization, and a recent development has 
been to combine these two effects in the context of a simple three-sphere model for 
the swimming of Chlamydomonas 50-52 . This simple model captures some of the 



important features of Chlamydomonas, namely, the ability to swim, the exact role of 
hydrodynamic interactions 50 , 52 , the existence of stable synchronized states and an 



emergent run-and-tumble behaviour which is observed when we add white noise to 
the driving force 51 . Here, we consider the model without added noise and explore 



its phase diagram and full parameter space to see when the model evolves into the 
synchronized state. We also investigate the stability of the synchronized state under 
various conditions, and the different types of behaviour that can be obtained form the 
model. These studies have revealed a number of intriguing features. 

2. The Model 



We arrange three beads in the x — y plane, each of radius a, on a frictionless scaffold, 
as shown in figure [T] We refer to the left, right and back beads with the subscripts '/', 
'r' and '6', respectively. Let Rq be the origin of the cell reference frame with respect 
to a lab frame. The cell axes x, y make an angle 9{t) with the lab axes X, Y. The left 
and right beads model the flagella and move on circular trajectories in the cell frame 



i = Ro - L^ + bhi, 


R/ 


r = Ro + Lx + bhr, 


R/ 


6 = R-o — Hy, 


Rfe 
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of radius b in opposite directions and with phases (pi and (pr] the back bead models the 
cell body and is fixed with respect to the cell frame. The positions and velocities of the 
beads are 

tio + Ley + b{^i-e)ii, (1) 

Ro - L9y + K0r. + 9)ir, (2) 

Ro - HOSc. (3) 

where the unit vectors fij and tj are in the normal and tangent directions of the circular 
trajectory of Rj. The left and right beads are driven by tangential forces Fj" and 
F^ respectively. Normal forces F" and F" are exerted by the beads in order to be 
constrained to the circular trajectories. The force on the back bead is such that the 
swimmer is force free and torque free: 

F, + F, + F, = 0, Ti + Tr + Tb = 0, (4) 

where Fj = F/tj + F"nj for i = l,r and Tj = Hj x Fj for j = b, /, r. 

The forces and velocities are related through hydrodynamic interactions between 
the beads: 

R, = 1f, + (G„-G„).F,-G„.F„ (5) 

Rr = -^Fr + {Grl — Grb) " F/ — Grb " F,., (6) 

Rft = -^{Fi + F,) + Gu ■ Fi + G,, ■ F„ (7) 

where ^ = Gnrja is the friction coefficient of each bead (77 is viscosity of the ambient fiuid). 

In the limit when a is small compared with all other length scales, the hydrodynamic 

interaction is described by the Oseen tensor Gjj = g^ |^ , [I+Tijfij) with r^j = Fj — r^- 

The phase difference 6 = (pr — 4>i evolves according to 

(Rr — Rb) ■ tj. — (Ri — R5) ■ t; + 

(8) 
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26 cos (6/2) {L coscj)- H sin < 



where = (0/ + 0r)/2 and 
• 1 



(Ri - Kb) ■ hi 



H cos (0 - 6/2) + L sin (0 - 6/2) 
(Rj. — Rfe) ■ rij. 



(9) 



H cos (0 + 6/2) + L sin (0 + 6/2) 

We solve equation [s] numerically for a choice of stroke pattern (driving forces) F^cpi), 
i = l,r. The long term behaviour of 6 depends on the stroke pattern and in many 
cases the initial condition. We compute 6, 9 and to leading order in a/L, since we 



do not require hydrodynamic interactions for synchronization 50 , but we include the 
next order hydrodynamic term when computing the velocities, since we need this second 
order affect to achieve a net swimming velocity. 
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3. Swimming velocity in the synchronized state 

First we consider the synchronized state where Fl{(j)) = F^{(j)) and 6 = 0, so that 
(pi = 4>r = (p- We do not worry about the stability of the synchronized state, 
which we consider in the next section, and assume that the swimmer stays in this 
state. Since the Reynolds number is low, we need to ask 'does the model achieve 
net propulsion?' If hydrodynamic interactions are not included, then the cell just 
moves forwards and backwards and there is no net motion. However, if we include 
hydrodynamic interactions, which vary in strength around the cycle, then the symmetry 
in the swimming stroke is broken and net propulsion is achieved. 

The magnitude and direction of the net swimming velocity depends on the ratios 
H/b and L/b. Figure [2] shows the net swimming velocity in the y direction for a range 
of {H/b, L/b) and constant driving force F^^cf)) = F^. For L > 2.25\H\, swimming is in 
the positive direction, otherwise the cell swims in the negative direction. Polotzek and 



Priedrich in reference 52 give the following explanation of why the cell may swim in 
either direction. The instantaneous velocity v = f/g is the ratio of the force — /y, which 
has to be applied to the back bead to prevent it from moving, and the friction coefficient 
g associated with towing the swimmer in the y direction. The force / oscillates during 
a stroke cycle and the hydrodynamic interactions, which reduce the magnitude of /, 
are strongest when the beads are closest together. On the other hand, the friction 
coefficient is largest when the beads are furthest apart and smallest when the beads are 
close together. The geometry of the swimmer determines which effect dominates and 
therefore whether the net swimming is in the positive or negative direction. Herein we 
fix the values a/L = 1/33, H/b = 1/5 and L/b = 2, which results in forwards swimming 
for Fo > 0. 

The the net velocity is also affected by the force profile and we consider driving 
forces F^{(f)) such that J^^ d(j)F^{(f)) = 2ttFo, where Fq is a fixed average force. The 
net velocity v can be written as v = 1/T j^ dtRb, where i?b = R;, ■ y, T = j^^ dcp/cp 
is the period, and we can write Jg dti?b = J^'' d(f)Rb{4>) / (p{4>) . In the synchronized 
state the force dependence cancels in the ratio Rh/cpi so the force only enters the net 
velocity expression through the 1/T term. In order to maximise the net velocity, we 
must minimise the period T = J^^ d(j) /</){(/)), where we can write 0(0) = -F*(0)<l>(0). 
Minimising T with the constraint J„'^d0F*(0) = 2ttFq tells us that a constant force 
profile F^{(f)) = Fq maximises the net velocity. Clearly, increasing Fq increases the net 
velocity. However as we shall see in the next section, the synchronized state is not stable 
when we choose a constant force profile. Friedrich et al. showed that a constant driving 
force can give a stable synchronized state if we change the direction of rotation of the 



beads, so Fq < 50 , 52 , which is equivalent to changing the sign of H, but here we 



choose to work with Fn > and H > 0. 
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Figure 2. Contour plot of dependence of velocity on the parameter lengths. The 
contours show the velocities in units of Fo/{67rr]b). The zero contour lies approximately 
on the line L = 2.25|iJ| for sufficiently large H/b. 



4. Synchronization and stability 

We consider force profiles of the forms F^( 



Fo(l + 6| sinn0i), where — 1 < a. 






= Foil 
< 1 and i = 



(n) 

a; cos n 



Lr. 



^ _, and Fl{(f)i) = 
The definition of the 



synchronized state that we use here is zero (or integer multiple of 27r) phase difference 
between the two fiagella, i.e. when S = 27m, n & Z. Initially we tried to analyst the the 
synchronization stability by linear stability analysis, but we are unable to do this because 
we cannot perform a valid Taylor expansion when = 0^, where H cos (ps + L sin 0^ = 0. 
We work numerically to avoid this Taylor expansion; for further details see the appendix. 
We identify five main types of stability of the synchronized state by looking at the 
evolution of 6 from different initial conditions (5(0o) = 0.1 for a number of 0o equally 
spaced in the range 0o ^ [0,27r): (i) All the initial conditions (5(0o) = 0.1 evolve into 
the synchronized state for all 0o (the synchronized state is stable), (ii) Some choices of 
00 evolve into the synchronized state and others choices evolve into an oscillating state, 
but there is a larger number of 0o that lead to synchronization than the number of 0o 
that leads to oscillations, (iii) Some choices of 0o evolve into the synchronized state and 
others choices evolve into an oscillating state; the numbers of 0o that lead to each type 
of behaviour are similar, (iv) Some choices of 0o evolve into the synchronized state and 
others choices evolve into an oscillating state, but there is a larger number of 0o that 
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(a) 



8 2tt 



(b) 




Figure 3. Evolution of S for driving force F*{(l}) = Fo(l + b^^^ smcj)) and initial 
condition 5{(j) = 0) = 1.3 with (a) &(i) = 0.7, (b) &(i) = 0.4, (c) b'^^'> = -0.1. The insets 
show a small part of the plot in more detail. The bottom axis (J)/2'k shows the number 
of cycles which increases monotonically with time. Bottom row: corresponding 0(0). 



lead to oscillations than the number of 0o that leads to synchronization, (v) All choices 
of 00 evolve into the oscillating state, (the synchronized state is unstable). 

Although the choice of initial condition 5q = 0.1 is arbitrary, we want to know 
how likely it is that a small perturbation from the synchronized state will decay back 
to synchronization, or whether it is likely to evolve into an oscillating state. This 
choice of 6o is suitable for this purpose. For type (v) stability, if we start in the 
synchronized state, then it is likely that some numerical noise will kick 6 into an 
oscillating state. It is possible for a small amount of noise to kick the synchronized 
state into an oscillating state for types (ii), (iii), (iv), with low probabihty for type (ii), 
then increasing probability for type (iii) and then type (iv). 



4.I. Equal beat patterns 

First we consider the case where Fl{(f)) = F^{(f)) = F^{(j)). For each choice of coefficient 
and initial condition, 6 evolves either to an integer multiple of 27t and remains at this 
value (synchronization); or it reaches a state where it oscillates about vr sinusoidally; or 
it reaches a periodic state near zero, but never reaching zero. Figure [3] shows examples 
of these three cases and the corresponding orientation 6. When 6 oscillates about tt, 
then the orientation oscillates about some fixed value. The cycle averaged motion is in 
a straight line, but the cell jiggles from side to side as well as backwards and forwards 
as it moves along. When b^^' = 0.4 and the oscillations are near zero, there is a net drift 
in the orientation so the net motion of the cell is along a curved trajectory. 

For many choices of -F*(0), the initial condition determines whether 6 evolves into 
the synchronized state or the oscillating state. Figure |4] shows the the dependence of S 
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Figure 4. Phase diagrams siiowing iiow S evolves for different initial conditions for 
F*(0) = 1 + a(i) cos with (a) a^^) = 0.4, (b) a'^' = 0.6, (c) a^^) = 0.8. Black squares 
represent initial conditions which lead to the synchronized state and white squares 
represent initial conditions which lead to a periodic oscillating state. The stability of 
the synchronized state is: (a) type (v) unstable; (b) type (i) stable; (c) type (ii). 



evolution on initial condition for -F*(0) = -fo(l + a^ ' cos(f)). A 63 x 69 grid is shown 
where each square represents an initial condition {(J)q,6o). A black square represents an 
initial condition for which 5 — )■ after a sufficiently long time; a white square represents 
an initial condition for which 6 continues to oscillate periodically as t — )■ oo. 

For a^^^ = 0.4, all initial conditions lead to an oscillating state and the synchronized 
state is unstable (type (v)). The black squares in figure |4] are initially in the synchronized 
state. Many squares along the line 6 = are white because a small amount of numerical 
noise drives the system away from the synchronized state. For a^^^ = 0.6, initial 
conditions close to 6 = lead to the synchronized state, but initial conditions far 
from 6 = lead to an oscillating state. The synchronized state is stable (type (i)). For 
a^^^ = 0.8, most initial conditions close to 5 = lead to the synchronized state, but 
a few initial conditions close to 5 = lead to an oscillating state and we have type 
(ii) stability; if 6 starts close to the synchronized state, it is likely to evolve into the 
synchronized state and it will stay in the synchronized state if there is no noise, but it 
is also possible for the cell to start close to the synchronized state and move away into 
an oscillating state. In an oscillating state the cell can still swim, but there will be more 
side to side movement. There appear to be a few white squares on the line 6o = 27r, 
however this is because the grid does not lie exactly on the 27r line, the grid contains 
points 6 = 6.2 and 6 = 6.3 and this small deviation from the synchronized state 6 = 27r 
is enough for evolution into an oscillating state for a few choices of (f)Q. 

Similar behaviour is observed for other harmonics, but with different ranges of 
coefficients giving the different stability types for the synchronized state. For example, 
figure |5] shows the 4th harmonic for three different coefficients of cosine and initial 
condition 6{(f) = 0) = 0.1. Figure [5[a) shows that 6 oscillate about vr for a^^^ = 0.4, 
oscillations are close to zero for a^^^ = 0.8 (it is interesting to note the 4 peaks in 
every cycle), and 6 evolves into the synchronized state for a^'^^ = 0.9. Figure k3I shows 
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Figure 5. Evolution of 6 for tlie 4tli liarmonic witli initial condition 6{4) = 0) = 0.1 
and coefficient (a) a^^) = 0.4, (b) a'") = 0.8, (c) a^^) = 0.9. (a), (b) S evolves into 
different oscillating states, (c) S evolves into the synchronized state. The insets show 
a small part of the plot in more detail. 
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Figure 6. Stabilities for range of coefficients and harmonics. Each horizontal band 
represents a harmonic, which is broken up into 49 blocks and each block represents 
a coefficient a*^"-'. Orange blocks represent stability (i); dark blue blocks represent 
stability (ii); pale blue blocks represent stability (iii); pale green blocks represent 
stability (iv); and bright green blocks represent stability (v). 



the stability of the synchronized state for the first 10 harmonics for a discrete range 
of coefficients when we choose the cosine term, i.e. -F/(</>i) = Fq{1 + a^"^ cos ncpi) . We 
see that there are more type (i) and (ii) force profiles for negative coefficients than for 
positive, showing that we are more likely to end up in the synchronized state when the 
coefficient is negative than when the coefficient is positive, for an arbitrary choice of 
harmonic. We see that for n > 2 there is no type (i) behaviour, so we cannot guarantee 
that we will reach the synchronized states for these higher harmonics. The unstable 
type (v) band around a^"'^ = gets narrower as n increases, so we only need weak phase 
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Figure 7. Initial condition phase diagrams for each type of stabiUty and force profile 
F*(0) — Fo{l + a*-"' cosn0). (a) Type (i) stable diagrams for n = 2 and a'") = —0.44. 
(b) Type (ii) stability for n = 7 and a*^"' = —0.8. (c) Type (iii) stability for n = 7 
and a(") = -0.44. (d) Type (iv) stability for n = 8 and a^") = -0.24. (e) Type (v) 
stability for n = 6 and a'"-* = —0.04. 



dependence to reach a synchronized state for higher harmonics, but the region of initial 
conditions which does lead to the synchronized state is very small. Even after reaching 
the synchronized state, it is likely that noise will drive S away from the synchronized 
state. Usually the stability moves towards the lower types of stability (more stable 
types) as \a^^^\ increases for each harmonic, although there are some exceptions. For 
example, when n = 1 we can start in a type (i) region, then as \a^'^^\ increases we move 
into a type (ii) region. We also see a type (iii) region surrounded by type (iv) regions 
on both sides for n = 3 and a^^^ > 0. 

For each type of stability shown in figure [6| we select an arbitrary force profile and 
show the full initial condition phase diagram in figure u\ In the type (i) stable case, we 
see that oscillating states can still evolve (see also figure 111), but the initial conditions 
for an oscillating state are not close to the synchronized case. In a few cases when we 
replace the cosine with sine, all initial conditions lead to the synchronized state, for 
example, the force profile i^*(0) = Fq(1 + b^^^sincp) for b^^^ > 0.6, except for a very 
thin dotted white curve through the middle of the phase diagram indicating an unstable 
oscillating state. However, if we look at this oscillating state for long enough, after some 
time the numerical noise causes it to move into the synchronized state. 
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We see from the bright green band down the center of figure |6] that we need some 
form of phase dependence in order to achieve synchronization. If we choose a constant 
driving force then S evolves into an oscillating state for all initial conditions. The 
synchronized state is unstable, so even if we start in the synchronized state, a small 
amount of numerical noise can drive the system into the oscillating state. Friedrich et 
al. showed that synchronization can occur with constant forcing when the direction 
of rotation of the beads is reversed (equivalent to if — ?■ —H and reversing swimming 
direction) but here we focus on the case F/ > and H > 0. 



The inspiration for our run-and-tumble model in reference 51 came from the 
initial condition phase diagrams. To see run-and-tumble we want to start in a stable 
synchronized state, then allow noise to move us temporarily into a white region of the 



phase diagram, before moving back into a black region. In reference 51 , we allowed 
noise to vary the coefficients a\ , so that a black square in the noiseless phase diagram 
can change to a white square when the instantaneous effect of the noise changes the 
value of a\ , then changes back to a black square when the instantaneous effect of the 
noise is smaller. This allows us to start in the synchronized state (with fluctuations due 
to the noise), then move away from the synchronized state when the instantaneous noise 
is large and we are at a suitable point in the phase diagram, then move back into the 



fluctuating synchronized state when the instantaneous noise is small. In reference 51 
we chose to work with the first harmonic and use a^^^ = 0.7 + C, where C is the noise 
term. We chose the value 0.7 because it is in the stability type (i) region, but lies close 
to the type (ii) region. In the fiuctuating synchronized state, when we move through 
the values of which are surrounded by white squares in the type (ii) phase diagram, 
there is the possibility to move into an oscillating state, but there are still plenty of 
black squares surrounding the line 5 = 0, so we can have long periods in the run phase. 
The noise causes fiuctuations of the position in the phase diagram, and this could 
also be a cause of run-and-tumble behaviour. For example, consider the phase diagram 
in figure [Wa). If we are in the synchronized state and noise is small enough such that 
fluctuations in S are within the black region, then tumbles will not occur. If the noise is 
larger, so 5 fluctuations move into the white region, then the cell could begin to move into 
the oscillating state, and after a few oscillations noise could kick the oscillations into the 
black region and the cell would move back towards a synchronized state. Elsewhere, we 
will consider the effects of adding noise to <p and 6, without any noise in the coefficient, 
to see if we obtain run-and-tumble behaviour this way and compare the statistics to the 
run-and-tumble obtained when noise is added to the coefficient. 

4-2. Mismatched coefficients 

If we start in the synchronized state, then 



dS 
d0 



5=0 \-^l ' ^r . 
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Figure 8. (a), (c), (e) Evolution of S and (b), (d), (f) corresponding orientation 
6 for force profiles Fl{4>i) = Fo{l + aiCOS(f>i) and initial conditions (5((/) = 0) = 0, 
6'(0 = 0) = and with (a), (b) -a; = Ur = 0.7; (c), (d) a/ = 0.7, a^ = 0.8; (e), (f) 
ai — —0.7, ar = 0.8. Each inset shows a small part of the main plot in more detail. 



which is non-zero for F/ ^ F^:, so the system does not stay in the synchronized 
state. We focus on the first harmonic and consider the case F/ = Fo(l + ai cos (pi) 
and F^: = Fq{1 + a^ cos0r) where o; ^ a^ (and we have dropped the upper index on the 
coefficient). When ai = —ar and \ai\ > 0.6, synchronization is frustrated and S oscillates 
about 2?™, n G 7j, shown in figure IsFa) for —ai = a^ = 0.7. The orientation of the cell 
drifts, shown in figure p[b) so the cell swims along a curved trajectory. For \ai\ < 0.6, 
the centre of 5 oscillations drifts away from the synchronized state, but remains close to 
27m. Swapping the signs of the coefficients swaps the direction of the orientation drift. 
When the coefficients have the same sign but different magnitudes, 5 oscillates about 
{2n + l)7r and the orientation oscillates about some fixed value. Figure [stc) shows the 
evolution of 5 for ai = 0.7, Oj. = 0.8. This type of behaviour is also seen for some 
choices of equal coefficients, for example, in figure p[e), (f) where b} = br = —0.1. 
Choices of coefficients which give this type of behaviour can be used to model a mutant 



of Chlamydomonas which swims with antiphase synchrony 54 . When the model swims 



with antiphase beating, the beat frequency is higher than when it swims with in phase 



beating, which has been observed in real Chlamydomonas cells 54 



When the coefficients have opposite signs and different magnitudes, there are two 
main types of behaviour that occur. For \ai\ < \ar\ and sgn(a;) = —1 or \ar\ < \ai\ and 
sgn(aj,) = —1, then 6 oscillates periodically and the corresponding orientation drifts 
in the negative direction in the former case and in the positive direction in the latter 
case. This is shown in figure |8](e), (f) for ai = — 0.7, a^ = 0.8. For \ai\ < \ar\ and 
sgn(a;) = +1 or \ar\ < \ai\ and sgn(ar) = +1, then oscillations in phase difference, 
S, drift in the negative (positive) direction and the orientation drifts in the positive 
(negative) direction in the former (latter) case. Figure M shows examples of this this 



Phase Dependent Forcing and Synchronization 



13 




(b) 
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Figure 9. Drifting S and corresponding orientation for (a), (b) F/ = Fo(l +O.2cos0;) 
andF^* = Fo(l-O.6cos0^); (c), (d) Fl = Fq {1-0 A cos (j)i) smd F^ = Fo(l + O.3cos0^). 
The insets show a small part of the plot in more detail. 



case where one bead completes more cycles than the other bead. There is a difference 
in the mean angular velocity of the two beads without choosing a different Fq for each 
bead. 



4-3. Combinations of harmonics 

So far we have considered force profiles with only one harmonic term. Now we 
consider driving forces with contributions from two harmonics. For simplicity we choose 
equal profiles for the left and right beads, -F/((/>) = F^{<P) = F'^{4>), and of the form 
F*{(j)) = -Fo(l + a^"^^ cosmcf) + a*-") cosn0), m ^ n and with the a*^*)'s chosen such that 
F*(0) > for all real 
2 and m 



m 



1, n 



Figure 10 shows the stability of the synchronized state for 
3. Each grid square represents a choice of driving 



l,n 



force with coefficients (a*^^\ a^-'-'), j = 2,3 and the colour represents the stability of the 
synchronized state for that particular driving force. We see that there are large regions 
for which the synchronized state is type (i) stable. 



5. Conclusion 



This simple mechanical model is able to evolve into a stable synchronized state for 
certain choices of parameters in the driving force when the initial condition is within 
some region of the synchronized state. We do not need hydrodynamic interactions to 
achieve stable synchronization; we include hydrodynamic friction on each bead, with 
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Figure 10. Phase diagrams sliowing stability of synchronized state for different 
combinations of parameters in driving force profiles (a) F*(0) = Fo{l + a^-^^ cos cj) + 
a^^) cos 20), (b) F*(0) = Fo(l + a^^' cos0 + a'^^^ cos 30). The colour of each square 
represents the following stability: Orange represents type (i); dark blue represents 
type (ii); pale blue represents type (iii); pale green represents type (iv); bright green 
represents type (v). 



force free and torque free conditions and a phase dependent driving force, which can be 
constructed with a suitable combination of harmonic terms. For many choices of force 
profile, some initial conditions allow the model to evolve into the synchronized state, 
while other initial conditions that are very close to the synchronized state lead to an 
oscillating state. There are some force profiles, including constant forcing, where there 
are no initial conditions that evolve into the synchronized state, and if the system starts 
in the synchronized state when the driving forces are equal, even a small amount of 
numerical noise can drive the system into an oscillating state. There are different types 
of periodic behaviour for different choices of parameters in the driving force; often the 
phase difference oscillates about vr, but sometimes the oscillations can occur close to 
zero with multiple peaks per cycle. 

When the parameter in the driving force is different for the left and right beads, we 
can get periodic oscillating states about a range of values, or we can get a drifting 
oscillating state, where one bead has a higher average angular velocity than the 
other. When the coefficients have equal magnitude and opposite sign, this can lead to 
oscillations about the synchronized state. This frustrated synchronization is interesting 
when we add intrinsic noise to the driving force, because then the behaviour of S is very 
similar for both opposite coefficients and equal coefficients, although the behaviour in 
the orientation is different in the two cases. 

The nonlinear mechanics of the system make it difficult to study analytically and 
it is not easy to predict the parameter ranges which give stable synchronization. Our 
numerical results have highlighted some of the main types of behaviour of the model. 

An important feature is that it is necessary to have some sort of phase dependent 
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driving force in order to have a stable type (i) or type (ii) synchronized state. When 
the phase dependence is only weak, then the synchronized state is unstable, which we 
see from figure [6] when the coefficient is small. The value of the coefficient at which 
the phase dependence becomes strong enough to give synchronization depends on the 
harmonic, whether we choose a positive or negative coefficient, and whether we choose 
sine or cosine. These latter choices are equivalent to adding a constant phase 0, n or 
±7r/2 in the harmonic term. 

This simple mechanical model shows a wide range of behaviour when we vary 
the parameter in the driving forces. The variety of stabilities suggests possibilities for 
developing run-and-tumble models, where noise can be used to to jump between regions 
of a phase diagram that lead to synchronization or oscillations, or jump between phase 



diagrams as we did in reference 51 
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Appendix 

We note that linear stability analysis can not be used to probe the stability of 
the synchronization, as we cannot perform a valid Taylor expansion when = 0^, 
where H cost's + LsiiKps = 0. For equal driving force profiles that we linearize, 
Fl{(j) - 6/2) ^ Fo{(f)) - Fi{(j))6/2, F^{(j) + 6/2) ^ Fo(0) + Fi(0)5/2, if we were to Taylor 
expand, then the linearized expression for 6' would be 

d5 ^ /(0;i^o,i^i) , . , . 

d0 ^ (i/cos0 + Lsin0)2 ' ^ ' ' 

which has a singularity at = 0s. The apparent singularity actually occurs at (pi = 0^ 
and at 0^ = 0s in the full expression, but the choice of constraining force ensures this 
zero in the denominator is canceled by the numerator. However, when we expand in 6 
and shift the singularity so that it occurs at = 0s, then the numerator is no longer 
zero at this point. The reason we have this zero in the denominator is the following: 
The torque free condition dll) is 

= (Ri - R,) X F, + (R, - Rb) X F„ 

= (F/(-Lcos0i + HsiiKpi + b) + F"(-Lsin0i - Hcoscpi) + 

F^{L cos (pr - iY sin0^ - b) + F;^{L siiKpr + //cos0,.))z. (A.2) 

Along with equations (1 3][5 7), we use (A.2) to solve for the constraining forces F", 



i = l,r. However, at 0j = 0s, -F" is multiplied by a term which vanishes, so the torque 
free condition can be satisfied without specifying FJ^. We over-constrain the system 
when we divide by zero and specify F" at 0j = 0s. Geometrically, 0s corresponds to 
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the phase where Rj — R;, is parallel to rij. Our numerical analysis of the full expression 
avoids this singularity. 
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